I. Soving the diffusion problem with explicit Euler integration & finite differences


In [337]:
using PyPlot
PyPlot.version


Out[337]:
v"1.3.1"

Problem description

For now, we want to solve the diffusion problem on a rectangular domain. The reaction part of the equation we want to solve in the end is currently set to zero. We therefore want to find a function $u$ defined as follows.

\begin{align} u: [a, b] \times [c, d] \times ℝ^+ &\to ℝ^+ \\ (x, y, t) &\to u(x, y, t) \end{align}

We therefore want to solve the following PDE:

$$ u_t = D Δ u = D (u_xx + u_yy) \\ u(a,y,t) = u(b,y,t) = u(x,c,t)=u(x,d,t)=0 $$

For this first solution, we use forward Euler integration in time and a second-order finite differences method in space. This leads to the following method to calculate the next time step:

$$ u^{n+1} = u^n + D \cdot \Delta t \cdot A u^n $$

Here, $A$ is the finite differences matrix for the 2D Laplace operator.

Define (rectangular) problem domain and resolution for solver


In [338]:
D = 0.01 # diffusion coefficient;

In [339]:
# define domain of the problem (rectangle) and time interval
xstart = 1
xend   = 3
ystart = 2
yend   = 5
tstart = 0
tend   = 10


Out[339]:
10

In [340]:
# define number of intervals in space
Nx = 64
Ny = 64


Out[340]:
64

In [341]:
dx = (xend - xstart) / Nx
dy = (yend - ystart) / Ny


Out[341]:
0.046875

In order to have a stable solution, we set $\Delta t$ based on $\Delta x$ and $\Delta y$, as follows:

$$ \Delta t \propto \frac{\Delta x^2 + \Delta y^2}{D} $$

where $D$ is the diffusion coefficient.


In [342]:
securityfactor = 10 # at least 9 is needed (with Nx=Ny=64)
Nt = int(securityfactor * (tend - tstart) * D / (dx^2 + dy^2))
dt = (tend - tstart) / Nt


Out[342]:
0.031746031746031744

Exact Solution

As a very simple scenario, we use the following initial conditions:

$$ u(x,y,t = 0) = C_0 \cdot sin(π \frac{x-a}{b-a}) \cdot sin(π \frac{y-c}{d-c}) $$

This leads to the following analytical solution:

$$ u(x,y,t) = C_0 \cdot sin( \frac{x-a}{b-a}) \cdot sin(π \frac{y-c}{d-c}) \cdot exp( -\left( \frac{1}{(b-a)^2} + \frac{1}{(d-c)^2} \right) \pi^2 D t) $$

In [343]:
C0 = 5
exact(x,y,t) = C0 * sin((x-xstart)/(xend-xstart) * pi) * sin((y-ystart)/(yend-ystart) * pi) *
    exp( - (1/(xend-xstart)^2 + 1/(yend-ystart)^2) * pi^2 * D * t)


Out[343]:
exact (generic function with 1 method)

In [344]:
C_exact = Float64[exact(x,y,tend) for x=xstart:dx:xend, y=ystart:dy:yend]
pcolormesh(C_exact)
colorbar()
clim((0,C0))


Create grid and set up initial conditions


In [345]:
# use exact function to set up initial conditions
C = Float64[exact(x,y,0) for x=xstart:dx:xend, y=ystart:dy:yend];

Define matrix for finite differences

We define the finite differences matrices seperately for the $\partial_{xx}$ and $\partial_{yy}$ operators. The vector of unknowns $u$ is organized in col-major fashion (neighboring elements have the same y-value but not the same x-value).

Thus, the matrix for $\partial_{xx}$ is block-diagonal and organized as follows:

\begin{align*} FD_x = \begin{pmatrix} B_1 & 0 & \cdots & 0 \\ 0 & B_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & 0 & B_{N_y} \end{pmatrix} & \hspace{1cm} \text{with} & B_i = \begin{pmatrix} -2 & 1 & 0 & \cdots & 0 \\ 1 & -2 & 1 & \ddots & \vdots \\ 0 & 1 & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & 1 \\ 0 & \cdots & 0 & 1 & -2 \end{pmatrix}_{N_x \times N_x} \end{align*}

The matrix for $\partial_{xx}$ has has non-zero entries in the diagonal and the $N_{xx}$-th upper and lower diagonals.

\begin{align*} FD_y = \begin{pmatrix} -2 & 0 & \cdots & 1 & \cdots & 0 \\ 0 & -2 & 0 & \ddots & \ddots & \vdots \\ \vdots & 0 & -2 & \ddots & \ddots & 1 \\ 1 & \ddots & \ddots & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & 0 \\ 0 & \cdots & 1 & \cdots & 0 & -2 \end{pmatrix}_{N_x N_y \times N_x N_y} \end{align*}

The sum of these matrices, combined with their step sizes, forms the finite difference matrix that is used for the integration.

$$ FD = \frac{1}{\Delta x^2} FD_x + \frac{1}{\Delta y^2} FD_y $$

In [346]:
# number of unknowns in x and y direction
ux = (Nx-1);
uy = (Ny-1);

In [347]:
# finite difference matrix for ∂u/∂x
FDx_local = spdiagm((ones(ux-1),-2*ones(ux),ones(ux-1)),(-1,0,1));
FDx = blkdiag({FDx_local for i in 1:uy}...);
# finite difference matrix for ∂u/∂y
FDy = spdiagm((ones(ux*(uy-1)),-2*ones(ux*uy),ones(ux*(uy-1))),(-ux,0,ux));

In [348]:
FD = 1/dx^2 * FDx + 1/dy^2 * FDy;

Do integration with forward Euler method


In [349]:
u = vcat(C[2:end-1,2:end-1]...);
iterationmatrix = dt * D * FD;

In [350]:
for i in 1:Nt
    u += iterationmatrix * u
end
C[2:end-1,2:end-1] = reshape(u,ux,uy)
pcolormesh(C)
colorbar()
clim((0,C0))


Compare to exact solution


In [351]:
vecnorm(C - C_exact)


Out[351]:
0.01457819261573872

In [352]:
# collected values for varying x (Ny = 64, Nt = 10000)
error_x = [
    0.3506375511215073    # Nx = 4
    0.12577467828850733   # Nx = 8
    0.045244709643395255  # Nx = 16
    0.016937758467202517  # Nx = 32
    0.0073057790525524555 # Nx = 64
    0.004444809104674145  # Nx = 128
]
Nxs = [4 8 16 32 64 128]'
loglog(Nxs, error_x)
grid(true)